Stable Vortex-Bright Soliton Structures in Two-Component Bose Einstein 

Condensates 



O 
< 

o 



CO 
Oh 

d 



> 

m 
oo 

o 
o 



K. J. H. Law, 1 P. G. Kevrekidis, 2 and Laurette S. Tuckerman 3 

1 Warwick Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK 
2 Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA 
3 PMMH-ESPCI, CNRS (UMR 7636), Univ. Paris 6 & 7, 75231 Paris Cedex 5, France 

We report the numerical realization of robust 2-component structures in 2d and 3d Bose-Einstein 
Condensates with non-trivial topological charge in one component. We identify a stable symbiotic 
state in which a higher-dimensional bright soliton exists even in a homogeneous setting with defo- 
cusing interactions, due to the effective potential created by a stable vortex in the other component. 
The resulting vortex-bright solitary waves, generalizations of the recently experimentally observed 
dark-bright solitons, are found to be very robust in both in the homogeneous medium and in the 
presence of parabolic and periodic external confinement. 

PACS numbers: 



Introduction. Vortices in nonlinear field theory have 
a time- honored history [![. They are among the most 
striking features of superfluids, play a role in critical cur- 
rent densities and resistances of type-II superconductors 
through their transport properties, and are associated 
with quantum turbulence in superfluid helium [2). The 
advent of Bose-Einstein condensates (BECs) 15 years ago 
0, nas produced an ideal setting for exploring rele- 
vant phenomena. Since the experimental observation of 
matter- wave vortices [5], by using a phase- imprinting 
method between two hyperfme spin states of a 87 Rb BEC 
[g[ , the road opened for an extensive examination of vor- 
tex formation, dynamics and interactions. Stirring the 
BECs [7] above a certain critical angular speed jil-[ll| 
led to the production of few vortices and even of 
very robust vortex lattices [l2|. These structures have 
been produced by other experimental techniques, such 
as dragging obstacles through the BEC [HI or the non- 




linear interference of condensate fragments [14|. Later, 
not only unit- c harg ed, but also higher-charged structures 
were produced [15[ and their dynamical (in) stability was 
examined. This field also has strong similarities and over- 
lap with the emergence of vortices and even vortex lat- 
tices in nonlinear optical settings; see e.g. lrl Il7|. 
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Another remarkable possibility in both BEC 



and in nonlinear optics, e.g. [2l|, is that of multi- 
component settings. Matter waves exhibit rich phase 
separation dynamics driven by the nonlinear interatomic 
interactions between different species or states that make 



up the BECs. Longitudinal spin waves [22] . transi- 
tions between triangular and interlace d squar e vortex lat- 
tices [23] , striated magnetic domains 
target patterns j^"" 



a square 



and robust 
have all been observed, as well as 



tunable interspecies interactions [27j and transitions be- 
tween miscible and immiscible dynamics (2^ . 

We interweave these two settings, motivated by [29| in 
which dark-bright solitons have been created in a quasi- 
ld two-component BECs. These structures were pre- 
dicted [loj and extended to more complex settings such 



FIG. 1: (Color online) The energetically stable S=l vortex- 
bright without external potential for R = 0.99, r ma x = 60 and 
N — 5900 (top row). The bottom row shows radial profiles 
of unit and higher charge vortex-bright solitons in a homoge- 
neous medium on a regular (left) and logarithmic (right) scale. 
All profiles are for N = 10000, R = 0.99, and r max = 80. 



as spinor condensates [31] (as dark-dark-bright, or dark- 
bright-bright solutions), but were only realized experi- 
mentally in 2008. These are often termed "symbiotic 
solitons", as the bright component would be impossi- 
ble to sustain under repulsive inter-atomic interactions 
(i.e., defocusing nonlinearities, as considered here), un- 
less the dark-component creates an "effective potential" , 
of which the bright soliton is a bound state. The coupled 
bright solitary waves [32] and the gap ones of [33] con- 
stitute additional examples of symbiotic structures. We 
consider higher-dimensional realizations [30[ i.e. vortex- 
bright solitons of various topological charge in 2d, as well 
as in 3d [49[ . We find these symbiotic configurations to be 
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robust, with or without parabolic external confinement. 
In an optical lattice, the unstable vortex may in fact be 
stabilized by the bright soliton. The stability persists in 
3d, while for traps elongated in the direction of the vortex 
core, additional negative energy (potentially instability 
bearing) modes [34| emerge, as in the single-component 
vortex [35j. The work of [5[ has already offered a proto- 
typical dynamical realization of such states (analogous to 
their quasi- Id counterparts of [30| by [29j| ) and attests 
to their experimental relevance. We will first give the 
physical setup, then discuss the numerical methods and 
lastly display the results, as well as future directions. 

Physical Setup. The non-dimensional Hamiltonian for 
a two-component condensate in the mean-field approxi- 
mation reads (36| : 

h = J ^(v^j^v^j+^Vw^+i^i 2 ^!^! 2 -^^^ 

(i) 

where \t(r) G C 2 is the pseudo-spinor order parameter, 
|*| 2 = (l^lM^t, M = diag{ M i,M2} is the diag- 
onal matrix of chemical potentials associated with the 
conservation of the number of atoms N± = J dr\^/i\ 2 
and N 2 = J dr\^ 2 \ 2 ; a related useful diagnostic is R = 
Ni/(Ni + 7V 2 ) = Ni/N. U is a 2 x 2 matrix accounting 
for the effectively nonlinear interatomic interactions. For 
the |1, — 1) and |2, 1) components of 8T Rb we can use [26| 
Un = 1.03, U 12 = U 2 i = 1 and U 22 = 0.97. These de- 
termine, through the negative sign of det(Z7) = \U\, the 
immiscible nature of the interactions leading to phase 
separation 19|, 26[. The confining potential is 



V(t,z) = ^|r| 2 + + A[sm 2 (2u r x) -f sin 2 (2cj r ?/)], 
4 4 v ' 



Vol 



(2) 

where Vmt is the parabolic component (often created 
magnetically) and ^ol the periodic (optical) lattice com- 
ponent. The time and length scales are l/w n and 
^h/muj n , where m is the atomic mass and uj n is an 
arbitrary frequency in Hz. For 87 Rb with scattering 
length of ai2 = 5.5nm, (uo ri uo z ) = 2tt x (8,40) Hz, 
and choosing uo n = 5/4cj 2 , the ratio between the ac- 
tual and non-dimensional number of atoms is N^ c ^d = 
(h/2muj n ) 3/2 (huj n / g 3d ) = 10, where g 3d = (4irh 2 ai 2 /m) 
is the dimensional interaction parameter. For a 2d reduc- 
tion, the interaction parameter is g 2 d = g3d(m(jO z /27rh) 1 ^ 2 
(e.g. [38]) and taking uj n = (j r , the amplification factor 
is 7Vf aC) 2d = 30. The equations of motion (\I/,c.c.) T = 
Ja(SH/S^,c.c.) T = JaDH, where J = diag(-iJ, ii") 
and a interchanges rows (3, 4) with (1,2), for this infinite- 
dimensional Hamiltonian system are 

W = -V 2 ^ + V(r)V + U\^\ 2 • * - MV, (3) 

The stability of stationary solutions is determined by 
the eigenvalues of the Hessian of the Hamiltonian, <jD 2 H, 



and of JcrD 2 H. Negative eigenvalues of <jD 2 H indi- 
cate energetic instability, since "dissipative" perturba- 
tions (e.g. from exchanges of atoms with the thermal 
cloud if the temperature deviates from zero) in the sys- 
tem can render them dynamically unstable, as can colli- 
sions with other eigendirections even in the pure Hamil- 
tonian (zero-temperature) system. The linear stability 
of the latter system is examined through the eigenvalues 
A = A r + i\i of JaD 2 H; instability arises when A r 7^ 
since, due to the Hamiltonian structure, the eigenvalues 
are symmetric over both the real and imaginary axes. 
From prior experience which is confirmed again here, lin- 
ear stability indicates evolutionary non-linear stability in 
the mean-field model, at least for time scales on the or- 
der of tens of seconds (this is not generically the case for 
non-linear static solutions of Hamiltonian systems). 

The variation can be posed in the or the 

{^reab ^imag} basis. The former is useful when the po- 
tential is axisymmetric, since then small excitations to a 
stationary solution \P = (^1 (r)e iSl ° , ^ 2 (r)e iS29 ) of the 
form ip = (ai,a 2 ) T (r)e xt + (61, 6 2 ) T (r)e A * £ will have def- 
inite angular momentum aj(r,0) = c\j(r)e lK(X i e . If we 



set K ai = k, then 



2Si, K a 



S\ + S 2 , and 



$i — $2, so a single index k will indicate the an- 
gular momentum of the excitation with given eigenvalue 
A. Hence, the spectrum of eigenvalues {A} can be decom- 
posed as the union of the spectra {X K } pertaining to an- 
gular momentum n. We will also assume 62 =0, so that 
Si = S and n a2 = n^. It has been shown numerically 
[341 ] and analytically [37[ that instability windows arise 
in a single component with topological charge S only for 
wave numbers with \k\ < S. The null eigenvalues cor- 
responding to gauge invariance appear in the spectrum 
of n = S. For a single component in a parabolic trap, 
an anomalous mode for n = S — 1 converges to zero as 
uo r — >• 0, accounting for translational invariance and lead- 
ing to the energetic stability of the S = 1 vortex without 
external potential. For each < k < S — 1 (5 > 1] an 
anomalous mode leads to windows of instability 0, E3] • 
We show that these can be significantly suppressed, al- 
though the 5—1 spectrum occasionally leads to small in- 
stability windows for a small fraction bright-soliton com- 
ponent, % < iVi, for large TV with parabolic trap (no 
windows were observed without the trap). 

Numerical Methods. Our methods extend those in, 
e.g., [39|, |4o[. The spatial discretization in (r, 0, z) em- 

ts Chebyshev polynomials to represent r dependence 
The Fourier modes representing and z make 
the Laplacian operators diagonal in these directions. To 
identify stationary states of (|3]), we first obtain an initial 
estimate via imaginary-time (i.e., replacing t — »• it) inte- 
gration using a first-order implicit /explicit Euler scheme 
with At = 10 -2 . We then refine the solution using New- 
ton's method. The linear system arising at each New- 
ton step is solved using the matrix-free IDR(s) algorithm 
which requires only the action of the Hessian. 



ton step 

0,111, 
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To accelerate inversion, we precondition the system with 
the inverse Laplacian, using its block diagonal struc- 
ture. Hence, we solve the system V -2 D 2 H (\I/ n ) A n = 
V~ 2 DH(^! n ) and update ^ n+i = ^ n - A n for n 

0. 1 Fewer than 5 Newton iterations usually achieve 

an accuracy of ||V" 2 Di^(^)|| Z 2/||^|| Z 2 < 1(T 12 . 



max [Re(X )],S=2 




FIG. 2: (Color online) Growth rate of the S — 1 mode as a 
function of R with N = 4000 for S = 1 (left), and growth rate 
of the S — 2 mode as a function of R x a±N for S = 2 (right). 

For each stationary solution we use the matrix- 
free Implicitly Restarted Arnoldi algorithm to iteratively 
compute the eigenpairs of the linearization JaD 2 H(ty) 
to a specified tolerance [47j. In order to find the de- 
sired eigenvalues we use inverse iteration, with the IDR(s) 
method and inverse Laplacian preconditioning to solve 
the linear systems, as above. Here, the preconditioner 
is taken to be [«/cr(V 2 )] _1 , so that each iteration solves 
V- 2 D 2 H{^)v n+1 = -V- 2 crJv n . 

We used a resolution in (r, z) of 40 x 64 x 80 to rep- 
resent non-axisymmetric solutions and eigenvectors. For 
axisymmetric solutions, quantitative accuracy requires 
only 30 radial modes for N < 1000, but up to 200 modes 
for larger N. For eigenvectors, we use only 5 + 1 modes 
in 6 (see introduction) and identify quantitatively all ex- 
pected invariant and negative directions and windows of 
instability from (34j . 

Results. In 2d (uo z — » oo) for uo r = A = 0, we first 
demonstrate in Fig. Q] the existence of an energetically 
stable (and hence also dynamically stable) vortex-bright 
soliton state. This is so for all of the R—N values that we 
have sampled. Notice the symbiotic nature of the state, 
as a bright soliton would be impossible to support un- 
der the repulsive/self-defocusing interactions considered 
herein. Indeed, a similar state exists for vortices of higher 
topological charge 5, as shown in Fig. [T]for 5 < 5. The 
logarithmic scale shows that the soliton is more localized 
for larger 5. In this case, the negative energy modes [34] 
in the spectra associated to < ft < 5 — 1 may lead to 
dynamical instability from complex quartets of eigenval- 
ues as a result of Hamiltonian-Hopf bifurcations which 
can occur upon collision with positive modes [44]. 

In a parabolic trap, the 5=1 vortex-bright structure 
remains dynamically stable. However, breaking of trans- 
lational invariance produces a negative energy mode in 
the k = 5 — 1 spectrum. Thus, an additional control 
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FIG. 3: (Color online) The evolution in time of the unsta- 
ble S = 2 solution for (R, N) = (0.93, 200) perturbed in the 
direction of the growing excitation. Top: approximate vortic- 
ity density iso-contours show trajectories of the two vortices. 
Bottom: first splitting and rejoining. 



parameter (through the atom number of the bright com- 
ponent) may lead to collisions of this mode with positive 
energy modes and, hence, rare isolated windows of in- 
stability arising for large R < 1 and large TV; see Fig. 
[2] (left) for such a window for 5 = 1 as a function of 
R when N = 4000. A similar feature has recently been 
shown in the dark-bright Id analog of the vortex-bright 
states [45| . 

For higher-charge vortex-bright structures, the same 
situation holds for the k = S — 1 spectrum, while for 
R = 1 windows of instability arise from the negative 
modes in the spectra of < k < S — 1. As R decreases, 
however, these windows of instability are generically sup- 
pressed by the increasing presence of the second bright 
component. Fig. [2] (right) depicts the growth rate of the 
5 — 2 spectrum for 5 = 2 over R—N parameter space. An 
example of an unstable solution perturbed in the growing 
excitation direction is depicted in Fig. [3l The vortices 
first split from the center, and begin to part and precess, 
but once they are far enough and the bright component is 
bimodal, they approach again and the bright component 
resumes uni- modality. The sequence repeats, similarly to 
single-component 5 = 2 vortices [481 ]. 

When we impose an additional sinusoidal lattice po- 
tential, A > 0, the one-component (R = 1) 5 = 1 vortex 
may become unstable (due to resonant eigenvalue colli- 
sions and ensuing oscillatory instabilities), at least for A 
sufficiently large [46]. The same holds for a large mass 
ratio R < 1. However, below a critical i?, once again the 
bright component has a stabilizing influence. 

The vortex-bright structure is stable in 3d without the 
trap and with periodic boundary conditions in z. Indeed 
this is immediately clear upon Fourier transforming in z, 
since the spectrum of the Hessian decouples into an infi- 




nite family of sub-spectra equal to the 2d spectra shifted 
by k 2 z , and hence it remains non- negative. It is stable in 
the trapped case as well for A = and uo z = 5cj r , and 
indeed for uo z > uo r . When uo z = uo ri the solution has an- 
other rotational invariance, and additional negative en- 
ergy modes emerge for uj z < uj r . For 2uj z = uo r there are 
at least two additional negative energy modes, although 
this may not lead to dynamical instability. Upon addi- 
tion of the lattice, A > 0,the results are expected to be 
similar to 2d (up to considerations of the aspect ratio of 
the harmonic trapping). See Fig. |4]for an example with 
uo z — 5cj r . 

Discussion. We have generalized the dark-bright, 
quasi- Id soliton that has been predicted theoretically and 
observed experimentally in BECs to that of a vortex- 
bright robust dynamical entity that emerges as a stable 
structure in both 2d and 3d condensates (although simi- 
lar concepts could be directly applicable to the nonlinear 
optics of defocusing optical media). We also examined 
relevant structures in the presence of parabolic (mag- 
netic) and periodic (optical) trapping and found that 
they remain stable. While instabilities may arise (e.g. for 
higher topological charge, S, or as a result of the lattice), 
these are usually alleviated/suppressed by the presence 
of the second component in 2d and 3d. 

It would be interesting to determine the robust exis- 
tence of such waveforms, which are well within the reach 
of recent experiments, e.g. 
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. Our study suggests 
that higher-charge vortices in a single component may be 
stabilized by an external blue-detuned laser-beam poten- 
tial acting as the bright-soliton here. Hence, the stabil- 
ity of such vortices should be systematically examined in 
the presence of external potentials. Other themes such 
as multi-vortex-bright soliton interactions and lattices 
would also be natural extensions of the present work. 
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